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ABSTRACT 



The evolution of supernova remnants (SNRs) is studied, with particular at- 
tention to the effect of magnetic fields with axisymmetric two-dimensional mag- 
netohydrodynamical simulations. The evolution of magnetic SNRs is the same 
as non-magnetic ones in the adiabatic Sedov stage. After a thin shell is formed, 
the shell is driven by the pressure of the hot interior gas (bubble). Evolution in 
the pressure-driven snow-plow phase is much affected by the magnetic field. The 
shell sweeping the magnetic field lines thickens owing to the magnetic pressure 
force. After 5 x 10 5 yr - 2 x 10 6 yr, the inner boundary of the thick shell begins 
to contract. This compresses the hot bubble radially and maintains its thermal 
pressure. Thus, the bubble forms a prolate spheroidal shape and becomes thin- 
ner and thinner, since it expands in a direction parallel to the magnetic field for 
Bq ^ 3/xG. Finally, the bubble contracts. The porosity of the hot low-density 
gas in ISM is reduced, taking the effect of the magnetic field into account. 
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nants 



1 hanayama@th. nao.ac.jp 

2 also School of Physical Sciences, The Graduate University for Advanced Studies (SOKENDAI); 
tomisaka@th.nao.ac.jp 



- 2 - 



1. Introduction 

Structure of the interstellar medium (ISM) is a key issue in the understanding of our 
Galaxy. A large fraction of the interstellar space is occupied with warm (T ~ 10 4 K) and 
hot (T ^ 3 x 10 5 K) media. The origin of the hot component is the energy released from 
supernovae and/or stellar winds. Shock waves of the supernova remnants (SNRs) hit the 
interstellar gas to heat it up and to sweep it into a dense shell. The evolution of SNRs 
has been studied for many years (for example, Cox 1972; Chevalier 1974; Cioffi et al. 1988). 
After the shock expansion speed decelerates below ~ 260n a08 kms~ 1 , the gas passed through 
the shock front is quickly cooled and a dense shell forms. The expansion of the shell is 
decelerated as the interstellar gas accumulates. Finally, the expansion speed of the shock 
front reaches the random speed of turbulent motion in the ISM and the shock wave decays 
to a sound wave. 

The interstellar magnetic field is seen at a glance to run in a direction parallel to 
the galactic disk. As for the strength, combining rotation and dispersion measures, Rand 
& Kulkarni (1989) obtained a turbulent field strength of ~ 5/iG and its cell size of ~55 
pc. The energy density of the magnetic field Bq/8ti ~ lCT 12 dyn cm~ 2 (.E?o/5pG) 2 attains 
that of the ISM turbulent motion v t with a density p as Po^ 2 /2 ^ lCT 12 dyn cur 2 (p /2 x 
10~ 24 gcm~ 3 )(t>i/10kms -1 ) 2 and an interstellar thermal energy p = kTp/fiH ~ 3 x 1CT 12 
dyncm _2 (T/10 4 K)(p /2 x lCT 24 gcm~ 3 ). These three energies are in approximate equiparti- 
tion. This seems to indicate the magnetic field plays an important role in the ISM. At the 
same time, the late evolutionary phase of SNRs must also be significantly affected by the 
interstellar magnetic field. 

Magnetic SNRs were first studied by Chevalier (1974) using spherically symmetric one- 
dimensional hydro simulations including only the magnetic pressure force. Slavin & Cox 
(1992) pointed out the effects of the magnetic field as (1) an SNR forms a thick shell sup- 
ported by the magnetic pressure force; (2) the hot gas heated by the shock wave in the 
adiabatic phase occupies the interior of the SNR. This hot cavity is effectively re-compressed 
in the later phase by the effect of the magnetic pressure from their spherically symmetric cal- 
culations. Similar one-dimensional simulations have also been done by Smith & Cox (2001). 
The tension force of the magnetic field has similar effects, broadening the shell and com- 
pressing the hot bubble, but is not included in these spherically symmetric simulations. At 
least two-dimensional magnetohydrodynamical (MHD) simulations are necessary. Balsara 
et al. (2001) have studied the evolution of SNRs expanding in the magnetic turbulent ISM 
using three-dimensional MHD simulations and they have explained the evolution of SNRs 
in the thermal X-ray and non-thermal radio continuum. However, their study is restricted 
to the adiabatic stage. In the present paper, we explore the total evolution of SNRs in the 
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magnetized ISM using axisymmetric MHD simulations. 

The late phase evolution of SNRs is essential to understanding of the structure of the 
ISM (Koo & Kang 2004). For example, hot ionized gas is believed to be confined in the 
cavities of SNRs, especially in large-volume old SNRs. However, the volumes of the hot 
cavities must be strongly affected by the interstellar magnetic field strength. We devote 
particular attention to the late-phase evolution of SNRs. 

The plan of the paper is as follows: in §2, model and numerical method are described. 
Section 3 shows numerical result. We pay special attention to the comparison between 
models with strong and weak magnetic field. It is shown that the evolution of SNRs is also 
affected by the ISM density and local heating rates. Expansion laws of the SNR shell and hot 
cavity are also shown in §3. Section 4 is devoted to discussion on how the volume fraction 
of the hot ISM is affected by the magnetic field strength. The reason why the shell density 
of SNR is determined is also discussed in 84. 



2. Model and Numerical Method 

We assume that a uniform interstellar gas with density no(cm~ 3 ) and magnetic field 
strength B (fiG) is extending. The half-width at half- maximum of the thick Gaussian com- 
ponent and the scale-height of the exponential component of the HI disk are equal to 265 
pc and 403 pc, respectively (Dickey & Lockman 1990). Since the size of the SNR is smaller 
than the scale-height of the HI, if H i, in this paper, we ignore the effect of the density strat- 
ification in the z-direction, which becomes important after the size of the bubble exceeds 
several times H m (Tomisaka 1998). The scale-height of the magnetic field is believed to be 
much larger than that of the gas. We assume here a uniform density and a uniform magnetic 
field strength. 

The basic equations are the ideal MHD equations. Using the cylindrical coordinates 
(z,r,(f>), assuming that we have no toroidal components of velocity and magnetic field = 
B<f, = and that variables are independent from the azimuth angle as d/d(p = 0, the basic 
equations are as follows: 
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where the variables have their ordinary meanings as p (density), v = (v z ,v r ) (velocity), 
p (pressure), T (temperature), B = (B z ,B r ) (magnetic flux density), and e = p/(7 — 
1) + pv 2 /2 (internal plus kinetic energy per unit volume). Equation (1) is the continuity 
equation; equations (2) and (3) are the equations of motion; equation (4) represents the 
energy conservation equation for internal plus kinetic energy. The induction equations for 
the poloidal magnetic fields are equations (5) and (6). The terms TqUh and A(T)nn 2 in 
equation (4) represent, respectively, the heating and radiative cooling rates, where % means 
the number density of the Hydrogen atoms. We use the radiative cooling rate estimated for 
the gas in the collisional ionization equilibrium with the solar metallicity (Raymond, Cox, 
& Smith 1976). 



2.1. Numerical Method 

We employ here the modified Lax-Wendroff scheme (Rubin & Burstein 1967), an explicit 
finite-difference method, to solve the basic equations (l)-(6). The number of grids and grid 
spacings are chosen as (N z , N r ) = (2500, 2500) and Az = Ar = 0.1 pc. Thus, our numerical 
domain covers a region of 250 pc x 250 pc. Artificial viscosity is added to smooth the jump 
at the shock front. For all the quantities U = (p, pv, e + B 2 /8n, B), we added an extra 
term to express artificial diffusion and solve equations like 

dU 

= V(cVC/), (7) 

where we take 
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with A = Az = Ar and A v = 4. The numerical scheme was tested by comparing known 
solutions which have been analytically and numerically obtained. We checked the program 
by comparing (1) the adiabatic SNR with the Sedov solution (Sedov 1959), (2) the spherically 
symmetric SNR with previously calculated (Slavin & Cox 1992) and (3) MHD shock tube 
problems (Brio & Wu 1988). The analytical solution of (1) is reproduced within a 5% relative 
error, for example, pressure and density distributions. The expansion law of the shock front 
of hydrodynamical SNRs of (2) coincides with previous work (such as that by Slavin & Cox 
(1992)) and the difference is not more than 5% for t <; lMyr. 

We assume a uniform ISM with number density no and temperature To. A uniform 
magnetic field is running in the z-direction as -Bo = (B z , B r ) = (Bq, 0). As for the radiative 
cooling rate, we use a spline function fit of cooling curves obtained by Raymond, Cox, & 
Smith (1976) (T > 10 4 K) and by Dalgarno & McCray (1972) (T < 10 4 K). We summarize 
the parameters used in Table 1. The heating coefficient T is chosen as follows: in models 
A-H (see Table 1), the heating rate I^nn is balanced by the cooling rate A(To)nn 2 in the 
ambient ISM. On the other hand, it is set to be 5 x 10 _27 ergs _1 in models AW-HW (see Table 
1) as minimum heating rate given by cosmic-ray (Spitzer 1978). We begin the simulation 
by adding thermal energy of E = 5 x 10 50 erg within the sphere of 2pc in radius. We 
assume a 10:1 H to He number ratio, which leads to a mean molecular weight for nuclei of 
u = 14/11 = 1.27 and that for fully-ionized ions and electrons of fi t = 14/23 = 0.61, and 
the number density of the Hydrogen atoms of % = (10/ll)n, where n means the number 
density of the gas. 



3. Results 

3.1. Model with B = 5uG 

In model A, we studied a magnetized ISM with B = 5/zG, n = 0.2cm~ 3 , and 
To = 10 4 K, which mimics the warm ionized ISM. With a magnetic pressure of Bq/8tt = 
10~ 12 dyncm~ 2 (_B /5/zG) 2 = 7.2 x 10 3 K cm~ 3 (_B /5/zG) 2 and a thermal pressure of p — 
(u /ut)n A;To = 5.8xl0~ 13 dyncm" 2 (n /0.2cm- 3 ) (T /10 4 K) = 4.2xl0 3 Kcm" 3 (n /0.2cm~ 3 ) 
(T /10 4 K), the plasma beta of the ISM is equal to f3 = 0.58(n /0.2cm- 3 ) (T /10 4 K) 
(i?o/5uG)~ 2 . In Figure 1, we plotted cross-cut views along the z- (a and c) and the r- 
axes (b and d). Upper panels (a and b) display density distributions, while the lower ones 
(c and d) display temperature distributions. The shell formation time at which the SNR 
changes its structure from the adiabatic Sedov phase (Sedov 1959) to the pressure-driven 
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snowplow phase is given (Cioffi et al. 1988) as 

The first snapshot (t = 7.97 x 10 4 yr; solid line) corresponds to the final stage of the adiabatic 
phase or the transition stage to the pressure- driven snow-plow phase. Figure 1 shows that 
the differences in the z- and r-cross-cuts are slight at this stage. A dip in the temperature 
distribution just after the shock front indicates that the radiative cooling begins to play a 
part at this stage. 

The second snapshot at t — 1 x 10 5 yr shows clearly that a thin shell is formed [Fig. 1(a) 
and (b)]. At this stage, Figure 1(a) and (b) indicates that the shell is thin and its widths 
propagating in z- and r-directions are similar. The temperature distributions also look 
similar to each other [Fig. 1(c) and (d)]. Peak density in the shell reaches maximum at this 
stage and declines gradually for the shell near the r-axis [Fig. 1(b)], while it increases further 
for the shell propagating near the z-axis [Fig. 1(a)] until t = 2.52 x 10 5 yr. 

At the age of 3.99 x 10 5 yr (dashed line), the SNR shows a shape far from spherical 
symmetry, that is, the shock front reaches R ~ 62 pc while Z ~ 57 pc. The shell width in 
the r-direction AR becomes much wider than that in the z-direction AZ (short-dashed line) 
as AR ~ 14 pc while AZ ~ 2 pc. The shape of the cavity which contains a hot low-density 
gas indicates an elongated shape along the z-axis. That is, the bubble (T ^ 10 5 K) half 
radius is equal to Z huh ~ 55pc while -Rtmb — 47pc. This difference in z- and r-directions 
comes from the effect of magnetic pressure which works to confine the gas in the r-direction. 
The post-shock thermal pressure (not shown) reaches ~ 4 x 10 4 K cm~ 3 on the z-axis, while it 
is equal to ~ 2 x 10 4 Kcm~ 3 for the shock propagating perpendicularly to the magnetic field 
(along the r-axis). In a case in which the magnetic field is parallel to the shock front, which 
is applicable to the shock expanding in the r-direction, the post-shock to pre-shock pressure 
ratio is given using the plasma beta of the pre-shock gas /3q and the Mach number of shock 
propagation speed Ai [the Rankin-Hugoniot relation, see eq.(A8)]. Since part of the pressure 
is attributed to the magnetic one, the post-shock thermal pressure of the magnetic shock 
must be lower than that of the non-magnetic one (the ratio of post-shock thermal pressure of 
the magnetic shock to that of the non-magnetic shock is plotted in Fig. 11). With the plasma 
P of the ISM p = 0.58 (model A) and the Mach number of the shock speed V s ~ 50km s _1 
of M. = V s /c s = 3.3 at t — 3.99 x 10 5 yr (c s denotes the adiabatic thermal speed), equation 
(A8) anticipates that the thermal pressure of magnetized post-shock gas is 0.36 times smaller 
than that of the non-magnetized one. This demonstrates that the post-shock magnetic field 
plays an essential role in determining the shock structure at this stage. 

The shell width continues to increase with time in both z- and r-directions. Even in 
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the z-direction, the shell width increases. This is understood as follows: In the case of 
hydrodynamic shock, the shell is confined under the effect of the post-shock pressure and 
the thermal pressure of the hot interior gas. Ram pressure working at the front of the shell, 
which is proportional to the shock speed squared, decreases with time, since the expansion 
of the shock front is decelerated. At the same time, the thermal pressure in the hot interior 
gas also decreases owing to the adiabatic expansion. These thicken the shell thick (Chevalier 
1974). In the r-direction, the magnetic force prevents the shell from being compressed and 
works to support the shell. The width of the shell increases owing to the decrease in the 
confining pressures and the difference between the shells propagating in z- and r-directions 
comes from the magnetic force . 

In Figures 2 and 3, we plotted density distributions, total (thermal + magnetic) pressure 
distributions and magnetic field lines. For each snapshot, the cross-cut views are shown in 
Figure 1. Figures 2(a) and 3(a) show the structure of SNRs at the final stage of the adiabatic 
phase (t = 7.97 x 10 4 yr). The interior pressure (mainly thermal) predominates over the ISM 
pressure and the SNR holds the spherical symmetry. The shell thickness is approximately 
equal to 10% of the radius, which decreases to form a cooled shell [Figs. 2(b) and 3(b); 
t — 1 x 10 5 yr]. The SNR looks spherically symmetric in the adiabatic phase, even with the 
magnetic field. Compared with panels (c) of t — 3.99 x 10 5 yr, it is shown that the shell width 
increases with time in the pressure- driven snow-plow phase after the cool shell formation, 
especially in the r-direction. 

As is shown in Figure 1, the distance of the shock front from the explosion site reaches 
R ~ 62pc in the r-direction, while Z ~ 57pc in the z-direction (panels c of t — 3.99 x 10 5 yr). 
This means that the shock wave propagates faster in the direction perpendicular to the 
magnetic field than in the parallel direction. Since the magnetic field lines are bent reaching 
the shock front after passing the shock front, the shock running outward is a fast shock. The 
phase velocity of the fast wave increases in the direction perpendicular to the magnetic field. 
It is natural that the fast shock propagates faster in the r-direction than in the z-direction. 
Thermal pressure in the hot low-density cavity is twice as high as the total pressure p+B 2 /8% 
of the ISM at this phase. However, it becomes lower than the total pressure of the ISM after 
t ^ 5 x 10 5 yr. Around this equilibrium epoch, the hot cavity begins to contract in the 
r-direction, while in the z-direction, the inner boundary of the shell continues to expand. 

Figures 2(d) and 3(d) show the structure at 2.52 x 10 6 yr, at which the snapshot is also 
made in Figure 1 (dash-dotted line). As indicated by a magnetic field line running nearest 
to the z-axis in these figures, the inner boundary of the shell contracts radially between 
3.99 x 10 5 yr (R ~ 47 pc) and 2.52 x 10 6 yr (R ~ 28 pc), while the outward-facing shock front 
continues to expand from R ~ 62pc at t — 3.99 x 10 5 yr to R ~ 145pc at t — 2.52 x 10 6 yr. 
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The portion traced by the above nearest magnetic field line corresponds to the temperature 
jump between a gas in the shell T ^ 10 4 K and a hot gas in the cavity T ^ 10 5 K. As for the 
expansion in the z-direction, the temperature transition located near Z ~ 55 pc (3.99 x 10 5 yr) 
expands to Z ~ 95 pc (2.52 x 10 6 yr) shown in Figure 1(c), while the cavity contracts in the 
r-direction as the radial size of the cavity decreases from R ~ 47 pc at t = 3.99 x 10 5 yr to 
R ~ 28 pc at t = 2.52 x 10 6 yr shown in Figure 1(d). This indicates that the expansion of 
the hot gas is prevented by the effect of the magnetic field in the direction perpendicular to 
the magnetic field. As for the shock front, panels (d) show that the angle measured from 
the shock normal vector to the pre-shock magnetic field line is positive for r ^ 40 pc. In 
other words, the magnetic field component parallel to the shock front does not change its 
sign passing through the shock front. This is characteristic of the fast shock, while in the 
region near the z-axis (r <; 40 pc), the parallel components have different signs comparing 
the pre- and post-shock magnetic fields. Thus, this part of the shock is the intermediate 
shock (Jeffrey & Taniuti 1964). 

The structure at t = 5.64 x 10 6 yr is shown in Figures 1 (long dashed line), 2(e) and 
3(e). The expansion of the front continues. In the thick diffuse shell propagating in the re- 
direction, the thermal pressure has an excess of only ~ 10% over the ISM and the excess total 
pressure in the shell is at most ~ 20% of the ISM, even including the magnetic pressure. The 
shock is weak. The expansion in the r-direction is steadily decelerated till t = 5.64 x 10 6 yr. 
Although the deceleration along the z-axis continues for t <^ 10 6 yr, the expansion speed 
becomes constant V ~ 20km s^" 1 after t ^ 10 6 yr. The shell propagating along the z-axis is 
not decelerated for t ^ 10 6 yr. This may seem strange. However, this can be explained as 
follows: The mechanical equilibrium in the r-direction requires that the thermal pressure in 
the hot low-density cavity p C3iV , where the magnetic field is weak, is in pressure equilibrium 
with the total pressure of the ISM as p cav ~ po + Bq/8tt in the late phase of SNRs. This leads 
to a thermal pressure imbalance of p cav > po- This difference of -Bq/87t pushes the shell near 
the z-axis, where no magnetic force works. The equation of motion for the part expanding 
along the z-axis is written 

where a represents the column density of a = f pdz. Constant force per unit area drives the 
mass-accumulating shell to move with a constant velocity as a oc Z oc t. 
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3.2. Anisotropy Driven by the Magnetic Field 

To see the departure from the spherically symmetric SNR more quantitatively, we define 
the ellipticity of the SNR shell as 

Ash = Z sh / R sh , (11) 

where and _R s h are the distances from the explosion site along the z- and r-axes and the 
distance is measured for the point at which the density changes from that of the ISM in the 
cross-cuts by more than 1% (Fig.l). As for the bubble, its ellipticity is calculated as 

Arab = Rbub/Zbub, (12) 

where Zb u b and i?bub are the distances from the surface of the bubble to the explosion site. 
The surface of the bubble is defined using the temperature isosurface of T = 10 5 K and the 
iso-density surface with the same density as the ISM n as 

-Rbub = min(i? T=10 5 K , R n=no ), (13) 

Zbub = min(Z T=10 5 K , Z n=no ), (14) 

where the former temperature isosurface indicates the surface in the pressure-driven phase 
(t ^ 10 5 yr) while the latter density isosurface does that in the adiabatic stage. We found the 
minimum of these two gives the bubble surface appropriately. It should be noted that these 
two factors are chosen to be smaller than unity for ordinary magnetized SNRs (Z sh < R sh 
while i?bub < Z huh ). Figure 4 plots ellipticities Ah and Aub- After t ^ 10 5 yr, Ah begins to 
decrease and reaches ~ 0.8 after t ^ 2 x 10 6 yr. This shows that the maximum difference 
between i? s h and is 20% and the SNR is observed as elongated in the r-direction (an 
oblate spheroidal shape), if we take the shock front. Asymmetry in the bubble shape also 
appears t ^ 10 5 yr and grows steadily. The bubble shrinks in the r-direction to form a 
very thin cigar- like structure Aub = -Rbub/^bub 1. The hot cavity in the pressure-driven 
expansion phase forms a prolate spheroid. 

Figure 4 also plots the ratio of the bubble radius to the shell radius in respective direc- 
tions as 

otz = Z bub /Z sh , (15) 

and 

Ot r = Rbub/Rsh- (16) 

Since in the Sedov phase the shell size is approximately equal to one-twelfth of the shock 
radius AR ~ R/12 (Sakashita & Ikeuchi 1996), a z and a r keep constant ~ 0.9 in the 
adiabatic phase t ^ 7 x 10 4 yr. As a cool shell forms, the difference between the shock radius 
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and that of the hot cavity decreases. The quantities a z and a r reach ~ 0.95 at the age 
of t ~ 10 5 yr. The shell width near the z-axis (AZ = — ^bub) continues to be thin for 
t <^ 4 x 10 5 yr and increases gradually after i ^ 4 x 10 5 yr. On the other hand, the shell 
width in the r-direction (AR = R^ — -Rtmb) increases steadily after t ;> 10 5 yr (a r decreases 
monotonically) . Figure 4 shows us that asymmetry from the spherically symmetric expansion 
begins just after the effective radiative cooling forms a thin shell. The diffuse shell formation 
owing to the magnetic field is triggered just after this shell formation. 

3.3. Model with B = 1//G 

In this section, we describe model C with a weak magnetic field of B = 1/xG. Other 
parameters are the same as model A. Since the magnetic pressure of the ISM is equal to 
B$/8ir = 4 x lO^dyncm-^Bo/l/uG) 2 , the plasma (3 of the ISM is as large as (3 = 14.5. 
Before and at the thin shell formation stage [Fig. 5(a) and (b)], the SNR is spherically 
symmetric and looks very similar to model A (first two panels of Figs. 2 and 3). However, 
although the structure of the magnetized SNR in B = 5fiG (model A) is far from being the 
spherically symmetric in the pressure-driven expansion stage [last three panels of Figs. 2 
and 3], the SNR in a weak magnetic field (model C) remains its shape spherically symmetric 
[Fig. 5(c) and (d)]. For example, at t — 1.42 x 10 6 yr, the shell extends in the r-direction 
from 65 to 89 pc and in the z-direction from 70 to 85 pc. 

In Figure 4 (b), we plot the ratio of bubble-to-shock radii, and major-to-minor axis 
ratio of the shock front and the bubble for model C. Compared with Figure 4 (a) of model 
A, although the ratio of bubble-to-shock radii is as small as a r ~ 0.5 at the age of 10 6 yr 
for model A, it is a r ~ 0.8 for model C. This means that the thickness of the shell, which 
propagates in the direction perpendicular to the magnetic field (r-direction), is much reduced 
in model C. The bubble begins to contract after 1.42 x 10 6 yr, which is much longer than that 
of model A t ~ 4 x 10 5 yr. The difference comes from the effect of the magnetic force working 
to compress the bubble. In both models A and C, a r decreases with time and reaches 0.1 in 
t~4x 10 6 yr (model A) and t~7x 10 6 yr (model C). 

As for the same ratio for the parallel (z-) direction, both models have similar values as 
a z ~ 0.93 for model A and a z ~ 0.88 for model C at the age of 10 6 yr. However, although 
model A remains at a z ~ 0.85 for the final phase t ^ 5 x 10 6 yr, a z continues to decrease 
in model C and a z ~ 0.1 for t ~ 7 x 10 6 yr. This means the shell propagating along the 
z-axis is confined to being thin in model A, while the shell thickness increases with time 
in model C. In the last paragraph of section 3.1, we already mentioned that the expansion 
of the shell propagating in the z-direction is driven by the thermal pressure in the bubble, 
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which is sustained by the radial compression due to the magnetic force. This seems to have 
the effect of keeping the shell that is propagating in the z-direction thin. 

3.4. Effect of the Density 

In this section, the models with high ISM density n = lcm~ 3 are shown. Model E 
has the same parameters as model A except for the interstellar density. The heating rate is 
taken to be 5 times larger than model A to keep T = 10 4 K. 

Equation (9) gives the timescale to enter the pressure- driven snow-plow phase from 
the Sedov stage as t sf ~ 3.1 x 10 4 yr(_E /5x 10 50 erg) 3/14 (n /lcm~ 3 )~ 4/7 , which is three 
times shorter than the models with low density n = 0.2cm -3 . We plotted the snapshot 
of the adiabatic stage at t — 2.83 x 10 4 yr in Figure 6 (solid line). The next snapshot is at 
t = 5.03 x 10 4 yr (dotted line), at which the density of the shell propagating in the r-direction 
takes the maximum. This shows the thin shell phase of SNR. 

At t = 5.03 x 10 5 yr (dashed line), i?bub takes the maximum R huh = 29pc and the bubble 
begins to contract in the r-direction. This timescale is almost the same as that of model A 
(t — 4 x 10 5 yr), even if the adiabatic to pressure-driven transition timescale t s { is three times 
shorter than low-density model A. At this stage, the asymmetry between z- and r-directions 
has developed. That is, the shell width of the r-direction AR = i? s h — -Rtmb — 15pc and that 
of the z-direction AZ = Z^ — Z^, ~ 5pc differ. And as a consequence, the shell propagating 
in the z-direction is thin, while that propagating in the r-direction is thick. 

The shock wave propagates outwardly in both z- and r-directions. The next snapshot 
is at t — 2.52 x 10 6 yr (dash-dotted line). The bubble has almost been crushed -Rtmb ^ lOpc. 
In ~ 3 x 10 6 yr the bubble is completely crushed. The contraction time-scale is much shorter 
than that of the low-density models (^6x 10 6 yr for model A). 

3.5. Effect of the Heating Rate 

In this section, we compare models with the same parameters except for the heating 
rate. In models AW-DW, we take To = 5 x 10~ 27 ergs -1 , which is ~ 1/10 smaller than 
previous models A-D. In Figure 7, we plot snapshots at t — 2.52 x 10 6 yr of the temperature 
distribution of models A (a) and AW (b) with B = 5fiG and models C (c) and CW (d) with 
B = l/iG. First, it is shown that the structures are very similar for models A and AW, 
even if r differs ~ 10 times. Although the height of the bubble Z huh differs several %, the 
radius of the bubble i?bub differs only slightly. On the other hand, the sizes of the bubbles 
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of models C and CW are different, as the bubble size differs ~ 20% for both R huh and Z huh . 

To closely see model A, as for the structure along the r-axis, the shell (30pc < r < 
140pc) in model A has a temperature of T ~ 10 4 K, while that for model AW is as low 
as T ~ 10 3 K owing to inefficient heating. As a result, the thermal pressure of model A 
[p th ~ 4 x 10 3 cm~ 3 K] is much higher than that of model AW, and decreases from the post- 
shock value of p t h — 4 x 10 3 cm~ 3 K to p t ^ ~ 3 x 10 2 cm~ 3 K near the bubble boundary. 
Although the thermal pressure differs one order of magnitude, the expansion of the bubble 
is unchanged. This is owing to the fact that magnetic pressure is dominant in the shell of 
these models. 

As for the structure along the z-axis, the temperature of the shell of model AW (T <^ 
10 3 K) is much lower than that of model A (T ~ 6 x 10 3 K). However, the thermal pressures 
in the shell look similar, as follows: p t h — 3 x 10 3 cm~ 3 K (model A) and p t ^ ~ 2 x 10 3 cm~ 3 K 
(model AW). Further, those in the bubble are equal to p t h — 8 x 10 3 cm -3 K (model A) and 
Pth - 6 x 10 3 cm~ 3 K (model AW). Since the thermal pressure in the bubble is driving the 
shell in the z-direction, this similarity explains the fact that the expansion is very similar 
even if T changes 10 times. As a result, it is concluded that the evolution does not depend 
on the heating rate for models with strong magnetic fields of B ^ 3fxG. 

In Figure 8, we plot the thermal pressure (a)-(b) and temperature distributions (c)-(d) 
of models C and CW. Before and at the shell formation stage (solid and dotted lines) the 
two evolutions are the same, since the heating rate has no effect on evolution at this stage. 
After the shell formation (dashed, dash-dotted, and long dashed lines), these two models 
show different evolutions. At t — 2.52 x 10 6 yr, the bubble occupies R <^ 60pc in model C 
while R <; 70pc in model CW, measured from the radius with a steep temperature gradient. 
At t — 5.64 x 10 6 yr, these are R <^ 35pc in model C and R 60pc in model CW (see also 
Fig. 9). This shows that the bubble shrinks to a greater degree in models with strong heating 
rates. 

Structures in the shell are significantly different: the thermal pressure and the temper- 
ature in model CW are much lower than those in model C. The low heating rate leads the 
shell to be cool and under low pressure. High pressure in the shell hinders the expansion of 
the bubble in model C and thus the bubble size is reduced. The fact that the pressure in 
the shell controls the bubble radius means the contraction of the bubble is partly owing to 
the gas pressure in the shell. As a result, we conclude that for models with weak magnetic 
field Bq <; 1/iG, the shell thickness increases and bubble size decreases as the heating rate 
increases. 

Comparing models E and EW with n = lcm -3 , even in the models with B =5fxG : the 
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size of the bubble depends on the heating rate as i?bub=26.3 pc (model E), 29.9 pc (model 
EW) at t = 7.97 x 10 5 yr. The radius of the bubble differs ~ 14% for i? bub between models 
E and EW, and the total pressure along the r-axis in model E is twice as high as that in 
model EW at t = 7.97 x 10 5 yr. The difference seems to come from a large difference in Yq 
(ro=2.5 x 10~ 25 erg s _1 for model E and 5 x 10~ 27 erg s _1 for model EW). In high-density 
models, the effect of T exceeds that of the magnetic field, and the evolution of the bubble 
is different in models E and EW, in contrast to models A and AW. 

3.6. Expansion Law 

In Figure 9, we plot the expansion law of the shell Z sn and i? sn for respective directions 
of z- and r-axes as well as that of the bubble Z bu b an d -Rbub for various models. Figure 9(a) 
shows that the expansion of the shock front -R sn is affected by the magnetic strength, if no 
and E are chosen to be the same. That is, the expansion is promoted with the strength of 
the magnetic field. This is related to the fact that the phase speed of the magneto-sonic wave 
propagating in the direction perpendicular to the magnetic field increases with the magnetic 
strength. However, the expansion of -R sn is independent of the value of the heating rate. 

As for the shell expanding along the z-axis [Fig.9(b)], the models with B ^ 3/xG show 
expansions different from those of from the models with weaker magnetic fields B ^ 3uG. 
The expansion of B ^ 3uG is promoted by the strength of the magnetic field, while the 
expansion speed in Z s h for B <^ 3uG finally approaches the adiabatic thermal speed c s . It has 
been pointed out that the late-phase expansion of Z sh with constant speed comes from the 
excess thermal pressure inside the bubble, which pushes the shell in the ^-direction (§3.1). 
Since the excess thermal pressure is well approximated by the ambient magnetic pressure as 
Bq/8tt ~ 10~ 12 dyncm~ 12 (i?o/5/iG) 2 , this excess becomes greater than the ambient thermal 
pressure p ~ 5.8 x 10- 13 dyncm- 12 (n /0.2cm- 3 )(T /10 4 K) only for B £ 3.8uG. This seems 
to explain the numerical result that the difference in Z sh disappears for B <^ 3/iG. 

Figure 9(c) and (d) shows the bubble sizes measured along the r- and z-axes, respec- 
tively. Figure 9(c) shows the expansion of -Rbub- This shows that the magnetic field suppresses 
the expansion of the bubble radius -Rbub, especially in the late phase evolution. This is a 
natural consequence of the fact that the contraction of the bubble is driven by the magnetic 
pressure and tension forces. Besides this, there is a clear difference in the expansions of 
models with strong heating (models A-D) and weak heating (models AW-DW). Especially 
in models B(W) - D(W) with Bq <^ 3«G, models with weak heating have larger radii -Rbub 
compared with those with strong heating. Since -R sn is not affected with To, this means that 
the shell becomes thinner in models with lower heating rates, for weak magnetic field models 
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Bq 3/iG. The fact that the shell density is higher in the weakly heated models BW-DW 
than in the strongly heated B-D is explained by the fact that higher density is necessary to 
equilibrate heating and cooling in the shell, which will be discussed in §4.2. 

Figure 9(d) shows the expansion of Z huh . There are two distinctly different evolutions of 
the bubble expansion in the z-direction. In models with B ^ 3/^G, Zb u b [models A(W) and 
B(W)] continues to expand. On the other hand, in models with B ^ 1/iG, Zbub expands 
first but contracts after (1.5 — 3) x 10 6 yr [models C(W) and D(W)]. Although all models 
C, CW, D, and DW show contraction in the final stage, the models with ineffective heating 
(CW and DW) have larger bubble sizes in the z-direction than those with effective heating 
(models C and D), respectively. This indicates models C and D have thicker shells than 
models C W and DW, respectively. This is similar to the shells expanding in the r-direction. 
In models with strong magnetic field, Z huh is larger for stronger magnetic field. This may 
seem strange, since the magnetic field has no force working in the ^-direction near the z-axis. 
However, as already seen, the bubble is contracting in the r-direction owing to the magnetic 
tension force for t ^ 10 6 yr for models A(W) and B(W) [Fig.9(c)]. The thermal pressure 
in the hot interior, which mainly drives the shell near the z-axis, increases with increasing 
magnetic field strength. This explains the fact that the inner boundary of the shell Zbub 
expands faster for stronger magnetic fields in the models with a non-negligible magnetic field 
of B £ 3^G. 



To understand the structure of ISM, the volume fraction of the hot gas fh is a key issue. 
The fraction of volumes covered by the hot cavities of SNRs is estimated by a quantity called 
porosity, which is defined as 



where tsn is the supernova rate per volume as ~ 10~ 13 yr _1 pc~ 3 (McKee & Ostriker 1977) 
and V(t) represents the volume covered by a hot cavity. We assume all the SNRs end their 
lives when the cavities contract and V(t) = 0. After obtaining the porosity, fh is estimated 
as fh = 1 — exp(— Q). In Figure 10, we plot the porosity and the volume fraction of the 
hot gas against the magnetic field strength. This shows clearly that the porosity decreases 
with the magnetic field strength. Comparing models A-D (diamonds), Q = 0.6r_i 3 for non- 
magnetic ISM while Q = 0.22r_ 13 for B = 5fiG. Here r_ 13 represents the SN rate per 
volume normalized by the standard value, 10 _13 yr _1 pc~ 3 . The anticipated volume filling 
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Discussion 



4.1. Porosity of ISM 
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factor of hot gas fh = 0.45 for non-magnetic ISM (model D) is reduced to f h = 0.2 for 
Bq = 5/iG (model A). In models with ineffective heating rates (models AW-DW), the bubble 
is compressed more weakly by the thick shell than models A-D. The porosity of these models 
is greater than that of models A-D. In these models, the effect of the magnetic field in 
reducing porosity is more remarkable, that is, Q = 1.95r_i3 for Bq = (model DW) while 
Q = 0.28r_i 3 for B = 5/xG (model AW). This anticipates f h = 0.86 (model DW) and 
f h = 0.24 (model AW). 

Even in models with ISM density of uq = lcm~ 3 , the effect of the magnetic field in 
reducing porosity is evident. Increasing Bq from 1/zG to 5/xG, porosity decreases by a factor 
of 1.7 from Q = 0.026r_ 13 to Q = 0.015r_ 13 . In the case where the average density of the 
ISM is high, the volume fraction of hot ISM is restricted unless the supernova rate is much 
higher than that of our Galaxy. Consequently, it is shown that the volume fraction of the 
hot gas is restricted as fh ~ 0.2 by the effect of the magnetic field, even if the average ISM 
density is as low as n ^ 0.2cm -3 . 

Observationally, in the solar neighborhood, the "not strongly absorbing" HI gas contains 
large holes, perhaps ranging up to 400 pc diameter, and the large holes occupy 10-20% of the 
volume, fh ^ 0.1 — 0.2 (Heiles 1980). For M31, another Sb galaxy, the filling factor which is 
derived for the HI holes, which are found in the range from 7 to 16 kpc from the nucleus, is 
about 1%, fh — 0.01 (Brinks & Bajaja 1986). For M33, an Sc galaxy, the filling factor for HI 
holes, allowing a scaling factor of 2 for non-detected holes, is less than 40%, fh <^ 0.4 (Deul 
& den Hartog 1990). HI holes arise by OB associations with multiple supernova explosions 
that produce supercavities, and are also produced by the independent, random individual 
supernova explosions (Heiles 1987). In this section, we considered the latter case, and the 
porosity and the filling factor of models for B = 5/iG presumably looks consistent with the 
observational data. 



4.2. Thermal Equilibrium 

To understand how the shell width is determined, we consider the gas in thermal equi- 
librium. Balance between the radiative cooling A(T)n 2 and heating T Q n leads 

m = «i, (18) 

i p 

where k represents the Boltzmann constant. The cooling function below T <^ 10 4 K is an 
increasing function of temperature and we assume A = A (T/10 4 K) 9 , ignoring the depen- 
dence of the ionization fraction x. Equation (18) gives the equilibrium pressure proportional 
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to 

p oc TqT 1 ^ (19) 
oc T 1/q n {q - 1)/q . (20) 

If we assume the post-shock gas is cooled and is in a thermal equilibrium [eq.(18)] under 
the post-shock pressure, it is shown that the density decreases with the deceleration of the 
shock velocity as 

n q/(q-l) q/(.q-l)y2q/( q -l) 

P ' 4 v s /r)1 x 

n oc — t-t. — -r oc — -r , (21) 

pi/(H) p Q 1 /(9- 1 ) v ' 

where we used the relation between the shock velocity V s and the post-shock gas pressure 
p. Since the power q = 3, the shell density n decreases as the shell expands with time, since 
the shock velocity V s decreases with time. 

In this discussion we have ignored the magnetic field, which works to prevent compres- 
sion. In the thick shell sweeping the magnetic field, dynamics is controlled by the magnetic 
field as shown in model A. In the extreme case, the change in density is controlled under 
the magnetic field. Although equation (18) does not hold near the shock front, the temper- 
ature and the thermal pressure reach the equilibrium values expected from equation (18) 
near the shell-bubble boundary. We calculate the equilibrium density under the pressure p 
in the shell, p eq (r) along the r-axis. The difference from this equilibrium density is given as 
A = (p(r) — p eq (r))/ p eq (r). The volume average of A decreases with time as 0.91,0.67,0.14, 
and 0.07 at t = 1 x 10 5 , 3.99 x 10 5 , 2.52 x 10 6 , and 5.64 x 10 6 yr respectively in model A. This 
shows that, with time, the density of the shell approaches the state of equilibrium expected 
from equation (18). 



Numerical calculations were carried out by VPP5000 in the Astronomical Data Analysis 
Center, National Astronomical Observatory of Japan. This work was partially supported by 
Grants-in-Aid for Scientific Research from MEXT, Japan (14540233,16036206). The authors 
would like to thank the CANS (Coordinated Astronomical Numerical Software) project who 
provided us with the MHD simulation code. 



-17- 



A. Rankin-Hugoniot Relation 

In the case that the shock front is proceeding in the direction perpendicular to the 
magnetic field (x-direction), the mass flux conservation advected with fluid motion is written 

as 

PlVxl = P2V X 2, (Al) 

where subscripts 1 and 2 denote the quantities of pre-shock and post-shock, respectively, and 
velocities are measured from the shock front which is moving with V s . Other momentum, 
energy, and magnetic flux conservations are written as 

B 2 B 2 
Pi + + pA = P2 + + P2V 2 x2 , (A2) 

( l 2 7 V v xl B\ (I 2 7 \ v x2 B 2 

v*i [^Piv xl + —Pi j + = v x2 ^-p 2 v x2 + — V2 J + — , (A3) 

v xl B x = v x2 B 2 . (A4) 



Using 



x s P5 = M = ^, (A5) 

Pi v x2 B l 

(A6) 



Pi 



the above equations (A1)-(A4) give following relations (Priest 1982; Sakashita & Ikeuchi 
1996) 



2(2-7) 2 M 
—fi-x +7 



( 1 -l)M 2 + - + 2 
Pi 



x - 7 (7 + 1) M 2 = 0, (A7) 



» = ^'( 1 -i)-^ i + 1 " (A8) 

where Ai = V s j (7P1/ ^Pi) 1 ^ 2 and /3i = 8npi/B 2 represent the Mach number of the shock 
speed against the pre-shock sound speed and the plasma beta for the pre-shock medium, 
respectively. 

We plot the ratio of the normalized post-shock pressure y to that of non- magnetized 
shock y(0i — > 00) as II = y((3i)/y((3i — > 00) against the shock Mach number Ai in Figure 
11. Even if M is fixed at 10, IT increases with (3 X as 11= 0.35, 0.67, 0.88, 0.96, and 0.99 for 
/3i— 0.1, 0.3, 1.0, 3.0, and 10, respectively. This indicates that small Ai and j3\ (slow shock 
speed and strong magnetic field) induce the ratio II to decrease, or in other words, the effect 
of the magnetic field becomes evident. 
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Using the Sedov solution, the expansion of the shock front is given as 

R s (t) = 1.15 I f/\ (A9) 

for the radius and 

(F \ 1/5 

V s (t) = 0.46 ( t~ 3/5 , (A10) 

for the expansion speed. These give the following typical values for R s and V s as 

r - - 3 *H*^T (^y'XvhT ' (au) 

V; = ISlkms-f £ °„ V /5 (^_)- /5 f ' V'*. (A12) 
V5xl0 50 ergy V0.2cm- 3 / V 10 y r / 

Using the Mach number of the shock expansion speed M = V s /c s and the plasma beta of the 
pre-shock ISM, the post-shock gas pressure is obtained in relation to the pre-shock pressure. 
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Table 1: Model parameters. 



Model 


n 


T 


Bo 


r 




(cm" 3 ) 


(K) 


fcxG) 


(ergs" 1 ) 


A 


0.2 


10 4 


5 


5 x 10" 26 


B 


0.2 


10 4 


3 


5 x 10- 26 


C 


0.2 


10 4 


1 


5 x 10~ 26 


D 


0.2 


10 4 





5 x 10~ 26 


AW 


0.2 


10 4 


5 


5 x 10- 27 


BW 


0.2 


10 4 


3 


5 x 10- 27 


cw 


0.2 


10 4 


1 


5 x 10~ 27 


DW 


2 


10 4 


o 


5 x 10- 27 


E 


1 


10 4 


5 


2.5 x 10~ 25 


F 


1 


10 4 


3 


2.5 x 10" 25 


G 


1 


10 4 


1 


2.5 x 10~ 25 


H 


1 


10 4 





2.5 x 10" 25 


EW 


1 


10 4 


5 


5 x 10- 27 


FW 


1 


10 4 


3 


5 x 10- 27 


GW 


1 


10 4 


1 


5 x 10- 27 


HW 


1 


10 4 





5 x 10- 27 
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Fig. 1. — Density distributions along the z- (upper-left: a) and r-axes (upper-right:b) of SNRs 
(model A). Those of temperature along the z- and r-axes are also shown in c (lower-left) and 
d (lower-right), respectively. Five snapshots are taken at the ages of t = 7.97 x 10 4 yr (solid 
line), 1 x 10 5 yr (dotted line), 3.99 x 10 5 yr (dashed line), 2.52 x 10 6 yr (dash-dotted line), and 
5.64 x 10 6 yr (long dashed line). 
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Fig. 2. — Structures of the magnetized SNR (model A). Density distributions (false color) 
and magnetic field lines (white lines) are shown. Five snapshots are taken at the ages of 
7.97 x 10 4 yr (a), 1 x 10 5 yr (b), 3.99 x 10 5 yr (c), 2.52 x 10 6 yr (d), and 5.64 x 10 6 yr (e). Panels 
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Fig. 3. — Total pressure p + B 2 /8n distribution and the magnetic field lines for model A. 
Five snapshots are taken at the ages of 7.97 x 10 4 yr (a), 1 x 10 5 yr (b), 3.99 x 10 5 yr (c), 
2.52 x 10 6 yr (d), and 5.64 x 10 6 yr (e). Panels (a)-(c) cover the region of 150pc x 150pc, 
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Fig. 4. — The ratios of the bubble radii to the shell radii as a z = Z huh /Z sh and a r = -Rbub/-R s h- 
We also plot the time variations of asymmetric factors /3 s h = Z sh / R sh and /3bub = -Rbub/^bub- 
Shell radii and are calculated as the distance of the shell from the explosion site 
along z- and r-axes, respectively. Bubble radii Zbub and -Rbub are the sizes of the hot interior 
cavity measured along the z- and r-axes. Models A (a: B = 5fiG) and C (b: B = 1/xG) 
are plotted. 
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Fig. 5. — Structures of the SNR expanding in the weak magnetic field model C. Density 
distributions (false color) and magnetic field lines (white lines) are shown. Four snapshots 
are taken at the ages of 7.97 x 10 4 yr (a), 1.59 x 10 5 yr (b), 1.42 x 10 6 yr (c), and 2.52 x 10 6 yr 
(d). 
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Fig. 6. — The same as Fig.l but for SNR in n = 1cm 3 (model E). Structures of the SNR 
expanding in the dense ISM with n = lcm~ 3 . Five snapshots are taken at the ages of 
2.83 x 10 4 yr (solid line), 5.03 x 10 4 yr (dotted line), 5.03 x 10 5 yr (dashed line), 2.52 x 10 6 yr 
(dash-dotted line), and 5.64 x 10 6 yr (long dashed line). 
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Fig. 7. — Comparison between models with strong heating rates and inefficient heating rates. 
Temperature distributions are plotted for models A (panel a), AW (panel b), C (panel c), 
and CW (panel d) at the age of 2.52 x 10 6 yr. 
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Fig. 8. — Comparison between models C and CW. Thermal pressure distribution along r- 
axis of model C (a) and that of model CW (b) are plotted. Temperature distribution along 
r-axis of model C (c) and that of model CW (d) are plotted. Snapshots are taken at the ages 
of 7.97 x 10 4 yr (solid line), 1.59 x 10 5 yr (dotted line), 1.42 x 10 6 yr (dashed line), 2.52 x 10 6 yr 
(dash-dotted line), and 5.64 x 10 6 yr (long dashed line). 
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Fig. 9. — Expansion laws of the shell size, _R s h (a) and Z sh (b) as well as the size of the 
bubble, i?bub (c) and Z huh (d) are plotted. Models A-D and AW-DW are shown. The figure 
shows that the expansion of the shell is affected by the strength of the magnetic field B , 
especially in R sh . A stronger magnetic field induces faster expansion, while, the bubble size 
(especially i?bub) is affected by both the strength of the magnetic field and the heating rate. 
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B (nG) B (^G) 

Fig. 10. — In panel (a), porosity of the volume occupied with the hot gas is plotted against the 
strength of the magnetic field. We have assumed the supernova rate of r§N = 10 _13 yr _1 pc -3 . 
In panel (b), an anticipated volume-filling factor of the hot gas fh is shown. Diamonds 
(models A-D), squares (models E-H), triangles (models AW-DW), and crosses (models EW- 
HW) represent models with different n and r . 
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Fig. 11. — Rankin- Hugoniot relation for the MHD shock propagating perpendicular to the 
magnetic field. The post-shock gas pressure normalized by that attained in the non-magnetic 
shock (II = y(Pi)/y(Pi — > oo)) is plotted against the shock Mach number {M). Each curve 
corresponds to different plasma (3. 



